function [By,Bx,u,xxi]=rfvar(ydata,lags,xdata,breaks)
% function [By,Bx,u,xxi]=rfvar(ydata,lags,xdata,breaks)
% This algorithm goes for accuracy without worrying about memory requirements.
% ydata:   dependent variable data matrix
% xdata:   exogenous variable data matrix
% lags:    number of lags
% breaks:  rows in ydata and xdata after which there is a break.  This allows for
%          discontinuities in the data (e.g. war years) and for the possibility of
%          adding dummy observations to implement a prior.
[T,nvar]=size(ydata);
[T2,nx]=size(xdata);
% note that x must be same length as y, even though first part of x will not be used.
% This is so that the lags parameter can be changed without reshaping the xdata matrix.
if T2 ~= T, disp('Mismatch of x and y data lengths'),end
if nargin<4
   nbreaks=0;breaks=[];
else
   nbreaks=length(breaks);
end 
X=zeros(T,nvar,lags);
if nargin<4
   nbreaks=0;breaks=[];
else
   nbreaks=length(breaks);
end
breaks=[0;breaks;T];
for nb=1:nbreaks+1
   for i=1:lags
      X(lags+breaks(nb)+1:breaks(nb+1),:,i)=ydata(lags+breaks(nb)+1-i:breaks(nb+1)-i,:);
   end
end
smpl=[];
for nb=1:nbreaks+1
   smpl=[smpl;[breaks(nb)+lags+1:breaks(nb+1)]'];
end
X=[X(:,:) xdata(:,:)];
X=X(smpl,:,:);
y=ydata(smpl,:);
[vl,d,vr]=svd(X(:,:),0);
di=1../diag(d);
B=vl'*y;
B=(vr.*repmat(di',nvar*lags+nx,1))*B;
u=y-X(:,:)*B;
xxi=vr.*repmat(di',nvar*lags+nx,1);
xxi=xxi*xxi';
B=reshape(B,[nvar*lags+nx,nvar]); % rhs variables, equations
By=B(1:nvar*lags,:);
By=reshape(By,nvar,lags,nvar);% variables, lags, equations
By=permute(By,[3,1,2]); %equations, variables, lags to match impulsdt.m
Bx=B(nvar*lags+(1:nx),:)';
